典型医学设计实验GEO数据分析 (step-by-step) - 数据获取到标准化
GEO是当今最大、最全的公共基因数据资源库，包括基因的表达、突变、修饰等信息，涵盖几乎所有的疾病，且单个实验检测样品数目较多，是我们分析、学习的很好资源。
实验设计
原始文章对14个溃疡性结肠炎病人 (Ulcerative colitis, UC)和15个克罗恩病病人 (Crohn’s disease, CD)的发炎组织和未发炎组织活检采样，用Affy芯片检测基因表达谱，研究发炎组织和未发炎组织的基因表达差异。(Genome-wide Pathway Analysis Using Gene Expression Data of Colonic Mucosa in Patients with Inflammatory Bowel Disease. Inflamm Bowel Dis. 杂志影响因子不高，但是领域的专业杂志)
检测并安装依赖的包
a = rownames(installed.packages())
install_bioc <- c("Biobase","oligoClasses","ArrayExpress", "pd.hugene.1.0.st.v1",
"hugene10sttranscriptcluster.db", "oligo", "arrayQualityMetrics",
"limma", "topGO", "ReactomePA", "clusterProfiler", "gplots", "ggplot2",
"geneplotter", "RColorBrewer", "pheatmap", "dplyr","stringr","genefilter")
for(i in install_bioc) {if(! i %in% a) BiocManager::install(i, update=F)}
# General Bioconductor packages
library(Biobase)
library(oligoClasses)
# Annotation and data import packages
library(ArrayExpress)
library(pd.hugene.1.0.st.v1)
library(hugene10sttranscriptcluster.db)
# Quality control and pre-processing packages
library(oligo)
library(arrayQualityMetrics)
# Analysis and statistics packages
library(limma)
library(topGO)
library(ReactomePA)
library(clusterProfiler)
# Plotting and color options packages
library(gplots)
library(ggplot2)
library(geneplotter)
library(RColorBrewer)
library(pheatmap)
# Formatting/documentation packages
#library(rmarkdown)
#library(BiocStyle)
library(dplyr)
#library(tidyr)
# Helpers:
library(stringr)
library(matrixStats)
library(genefilter)
#library(openxlsx)
#library(devtools)从ArrayExpress下载原始数据
实验检测采用的是Affy芯片 (A-AFFY-141 - Affymetrix GeneChip Human Gene
1.0 ST Array [HuGene-1_0-st-v1])，原始数据为CEL格式，存储于ArrayExpress，索引号是E-MTAB-2967。
使用getAE函数下载原始数据和注释数据到当前工作目录
(Rmd所在目录)，返回一个包含所有下载文件名的列表。
# type: 'raw' to download and extract only the raw data, 'processed' to download and extract only the processed data or 'full' to have both raw and processed data.
# anno_AE <- getAE("E-MTAB-2967", type = "raw")若自动下载不成功，手动点击下载图中File部分所有文件，放置到当前目录。(也可后台回复affy数据 获取)
# type: 'raw' to download and extract only the raw data, 'processed' to download and extract only the processed data or 'full' to have both raw and processed data.
# 修改参数local=T，读入当前目录的数据
anno_AE <- getAE("E-MTAB-2967", type = "raw", local=T)
## Warning in getAE("E-MTAB-2967", type = "raw", local = T): No processed data
## files found in directory /disk1/train/GEO/f1000anno_AE里面存储了所有的文件的路径和名字信息。并把原始数据解压缩，释放出CEL文件到当前目录，方便后续读取。
anno_AE
## $path
## [1] "/disk1/train/GEO/f1000"
##
## $rawFiles
## [1] "164_I_.CEL" "164_II.CEL" "183_I.CEL" "183_II.CEL" "2114_I.CEL"
## [6] "2114_II.CEL" "2209_A.CEL" "2209_B.CEL" "2255_I.CEL" "2255_II.CEL"
## [11] "2400_I.CEL" "2400_II.CEL" "2424_A.CEL" "2424_B.CEL" "255_I.CEL"
## [16] "255_II.CEL" "2826_I.CEL" "2826_II.CEL" "2853_I.CEL" "2853_II.CEL"
## [21] "2978_I.CEL" "2978_II.CEL" "2987_I.CEL" "2987_II.CEL" "2992_I.CEL"
## [26] "2992_II.CEL" "2995_I.CEL" "2995_II.CEL" "321_I.CEL" "321_II.CEL"
## [31] "3222_I.CEL" "3222_II.CEL" "3223_I.CEL" "3223_II.CEL" "3226_I.CEL"
## [36] "3226_II.CEL" "3233_I.CEL" "3233_II.CEL" "3258_I.CEL" "3258_II.CEL"
## [41] "3259_I.CEL" "3259_II.CEL" "3262_I.CEL" "3262_II.CEL" "3266_I.CEL"
## [46] "3266_II.CEL" "3269_I.CEL" "3269_II.CEL" "3271_I.CEL" "3271_II.CEL"
## [51] "3302_I.CEL" "3302_II.CEL" "3332_I.CEL" "3332_II.CEL" "848_A.CEL"
## [56] "848_B.CEL" "888_I.CEL" "888_II.CEL"
##
## $rawArchive
## [1] "E-MTAB-2967.raw.1.zip" "E-MTAB-2967.raw.2.zip"
##
## $processedFiles
## NULL
##
## $processedArchive
## character(0)
##
## $sdrf
## [1] "E-MTAB-2967.sdrf.txt"
##
## $idf
## [1] "E-MTAB-2967.idf.txt"
##
## $adf
## [1] "A-AFFY-141.adf.txt"原始数据解释
ArrayExpress的每个数据根据MAGE-TAB (MicroArray Gene Expression
Tabular)指南存储，包含5种不同类型的文件：
IDF (Investigation Description Format)：研究描述文件，包含实验信息如题目、描述、提交者联系方式、实验操作过程等。
ADF (Array Design Format): 芯片设计文件
SDRF (Sample and Data Relationship Format): 样本属性文件，如实验分组、处理方式、取样部位等。
原始数据文件 (raw)
加工后的数据文件 (processed data files)
ExpressionSets数据结构描述
组学数据通常比较复杂，包含很多不同的部分，如实验样品信息、基因组注释信息和实验数据，对应到芯片数据是样品属性和分组信息、不同来源的基因ID信息和功能注释信息、基因表达矩阵。
为了更好的组织这些数据，Bioconductor的Biobase包定义了一个标准化的数据结构
(ExpressionSet类)存储这些数据。ExpressionSet类包含下面几部分：
assayData: 芯片实验的表达数据，探针在行，样品名字在列，用函数
exprs获取metaData
phenoData: 样品描述信息，样品名字在行，实验分组、处理方式、取样部位在列。通常是
SDRF文件的信息的读入。用函数pData获取。featureData: 基因注释特征信息，行为探针名字，列为探针对应的基因、转录本名字和相应的注释信息等。用函数
fData获取。experimentData: 对实验描述的其它信息。
ExpressionSet类可以帮我们协调数据修改、提取过程中的所有信息的一致性。但是要注意phenoData的行名字必须与assayData的列名字一致，都是样品标识符；assayData的行名字必须与featureData的行名字一致，都是基因标识符，具体见下图：
导入数据，存储为”ExpressionSet”
读入SDRF数据。
sdrf_location <- file.path("E-MTAB-2967.sdrf.txt")
SDRF <- read.delim(sdrf_location)我们查看下，读进来的数据长什么样子，有哪些信息？
SDRF[1:3,1:5]总共有哪些样品相关的信息？取样个体标记、物种、疾病、取样部分、是否患病、raw data存储的位置和名字等。
t(SDRF[1, ])
## 1
## Source.Name "164_I"
## Characteristics.individual. "164"
## Characteristics.organism. "Homo sapiens"
## Characteristics.disease. "Crohn's disease"
## Characteristics.organism.part. "colon"
## Characteristics.phenotype. "non-inflamed colonic mucosa"
## Material.Type "organism part"
## Protocol.REF "P-MTAB-41361"
## Protocol.REF.1 "P-MTAB-41363"
## Extract.Name "164_I"
## Protocol.REF.2 "P-MTAB-41364"
## Labeled.Extract.Name "164_I:Biotin"
## Label "biotin "
## Protocol.REF.3 "P-MTAB-41366"
## Assay.Name "164_I_"
## Technology.Type "array assay"
## Array.Design.REF "A-AFFY-141"
## Term.Source.REF "ArrayExpress"
## Protocol.REF.4 "P-MTAB-41367"
## Array.Data.File "164_I_.CEL"
## Comment..ArrayExpress.FTP.file. "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip"
## Factor.Value.disease. "Crohn's disease"
## Factor.Value.phenotype. "non-inflamed colonic mucosa"Array.Data.File存储了样品名字信息
(CEL文件名)，用作行名字。把SDRF数据表转为AnnotatedDataFrame格式用于下游构建ExpressionSet对象。
rownames(SDRF) <- SDRF$Array.Data.File
SDRF <- AnnotatedDataFrame(SDRF)Affymetrix芯片的原始数据是CEL文件，里面包含检测到的探针杂交密度值，代表原始基因表达量。另外每个CEL还含有额外信息，如芯片类型、扫描时间等，常用于做批次校正。
oligo包的read.celfiles函数可以读取这些文件。(虽然不影响，但有几个文件命名不规律，如164_I_多写了个下划线，有几个样品用A,B代替了I,II。不规范的名字是不太利于批量分析和下游识别的，引以为戒。)
raw_data <- oligo::read.celfiles(filenames = as.character(SDRF$Array.Data.File),
verbose = FALSE, phenoData = SDRF)
stopifnot(validObject(raw_data))获得的raw_data就是一个ExpressionSet对象。raw_data@phenoData@data(或者pData(raw_data))是对应的样品属性信息，与SDRF信息一致。raw_data@annotation获取注释平台信息pd.hugene.1.0.st.v1。
str(raw_data)
## Formal class 'GeneFeatureSet' [package "oligoClasses"] with 9 slots
## ..@ manufacturer : chr "Affymetrix"
## ..@ intensityFile : chr NA
## ..@ assayData :<environment: 0x815aba0>
## ..@ phenoData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 23 obs. of 2 variables:
## .. .. .. ..$ labelDescription: chr [1:23] NA NA NA NA ...
## .. .. .. ..$ channel : Factor w/ 2 levels "exprs","_ALL_": 2 2 2 2 2 2 2 2 2 2 ...
## .. .. ..@ data :'data.frame': 58 obs. of 23 variables:
## .. .. .. ..$ Source.Name : Factor w/ 58 levels "164_I","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Characteristics.individual. : int [1:58] 164 164 183 183 2114 2114 2209 2209 2255 2255 ...
## .. .. .. ..$ Characteristics.organism. : Factor w/ 1 level "Homo sapiens": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.disease. : Factor w/ 2 levels "Crohn's disease",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.organism.part. : Factor w/ 1 level "colon": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Characteristics.phenotype. : Factor w/ 2 levels "inflamed colonic mucosa",..: 2 1 2 1 2 1 2 1 2 1 ...
## .. .. .. ..$ Material.Type : Factor w/ 1 level "organism part": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF : Factor w/ 1 level "P-MTAB-41361": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.1 : Factor w/ 1 level "P-MTAB-41363": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Extract.Name : Factor w/ 58 levels "164_I","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Protocol.REF.2 : Factor w/ 1 level "P-MTAB-41364": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Labeled.Extract.Name : Factor w/ 58 levels "164_I:Biotin",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Label : Factor w/ 1 level "biotin ": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.3 : Factor w/ 1 level "P-MTAB-41366": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Assay.Name : Factor w/ 58 levels "164_I_","164_II",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Technology.Type : Factor w/ 1 level "array assay": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Array.Design.REF : Factor w/ 1 level "A-AFFY-141": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Term.Source.REF : Factor w/ 1 level "ArrayExpress": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Protocol.REF.4 : Factor w/ 1 level "P-MTAB-41367": 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Array.Data.File : Factor w/ 58 levels "164_I_.CEL","164_II.CEL",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ Comment..ArrayExpress.FTP.file.: Factor w/ 2 levels "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Factor.Value.disease. : Factor w/ 2 levels "Crohn's disease",..: 1 1 1 1 1 1 1 1 1 1 ...
## .. .. .. ..$ Factor.Value.phenotype. : Factor w/ 2 levels "inflamed colonic mucosa",..: 2 1 2 1 2 1 2 1 2 1 ...
## .. .. ..@ dimLabels : chr [1:2] "rowNames" "columnNames"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ featureData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 0 obs. of 1 variable:
## .. .. .. ..$ labelDescription: chr(0)
## .. .. ..@ data :'data.frame': 1102500 obs. of 0 variables
## .. .. ..@ dimLabels : chr [1:2] "featureNames" "featureColumns"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ experimentData :Formal class 'MIAME' [package "Biobase"] with 13 slots
## .. .. ..@ name : chr ""
## .. .. ..@ lab : chr ""
## .. .. ..@ contact : chr ""
## .. .. ..@ title : chr ""
## .. .. ..@ abstract : chr ""
## .. .. ..@ url : chr ""
## .. .. ..@ pubMedIds : chr ""
## .. .. ..@ samples : list()
## .. .. ..@ hybridizations : list()
## .. .. ..@ normControls : list()
## .. .. ..@ preprocessing : list()
## .. .. ..@ other : list()
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 2
## .. .. .. .. .. ..$ : int [1:3] 1 0 0
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ annotation : chr "pd.hugene.1.0.st.v1"
## ..@ protocolData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 2 obs. of 2 variables:
## .. .. .. ..$ labelDescription: chr [1:2] "Names of files used in 'exprs'" "Run dates for files used in 'exprs'"
## .. .. .. ..$ channel : Factor w/ 2 levels "exprs","_ALL_": 2 2
## .. .. ..@ data :'data.frame': 58 obs. of 2 variables:
## .. .. .. ..$ exprs: Factor w/ 58 levels "164_I_.CEL","164_II.CEL",..: 1 2 3 4 5 6 7 8 9 10 ...
## .. .. .. ..$ dates: Factor w/ 58 levels "2010-06-09T11:39:13Z",..: 56 55 9 10 6 8 58 57 5 7 ...
## .. .. ..@ dimLabels : chr [1:2] "rowNames" "columnNames"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. ..@ .Data:List of 4
## .. .. .. ..$ : int [1:3] 3 5 1
## .. .. .. ..$ : int [1:3] 2 42 0
## .. .. .. ..$ : int [1:3] 1 3 0
## .. .. .. ..$ : int [1:3] 1 0 0
# 如果用Rstudio, View(raw_data)显示更好
# View(raw_data)查看样本属性信息
head(Biobase::pData(raw_data))
## Source.Name Characteristics.individual. Characteristics.organism.
## 164_I_.CEL 164_I 164 Homo sapiens
## 164_II.CEL 164_II 164 Homo sapiens
## 183_I.CEL 183_I 183 Homo sapiens
## 183_II.CEL 183_II 183 Homo sapiens
## 2114_I.CEL 2114_I 2114 Homo sapiens
## 2114_II.CEL 2114_II 2114 Homo sapiens
## Characteristics.disease. Characteristics.organism.part.
## 164_I_.CEL Crohn's disease colon
## 164_II.CEL Crohn's disease colon
## 183_I.CEL Crohn's disease colon
## 183_II.CEL Crohn's disease colon
## 2114_I.CEL Crohn's disease colon
## 2114_II.CEL Crohn's disease colon
## Characteristics.phenotype. Material.Type Protocol.REF
## 164_I_.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 164_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## 183_I.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 183_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## 2114_I.CEL non-inflamed colonic mucosa organism part P-MTAB-41361
## 2114_II.CEL inflamed colonic mucosa organism part P-MTAB-41361
## Protocol.REF.1 Extract.Name Protocol.REF.2 Labeled.Extract.Name
## 164_I_.CEL P-MTAB-41363 164_I P-MTAB-41364 164_I:Biotin
## 164_II.CEL P-MTAB-41363 164_II P-MTAB-41364 164_II:Biotin
## 183_I.CEL P-MTAB-41363 183_I P-MTAB-41364 183_I:Biotin
## 183_II.CEL P-MTAB-41363 183_II P-MTAB-41364 183_II:Biotin
## 2114_I.CEL P-MTAB-41363 2114_I P-MTAB-41364 2114_I:Biotin
## 2114_II.CEL P-MTAB-41363 2114_II P-MTAB-41364 2114_II:Biotin
## Label Protocol.REF.3 Assay.Name Technology.Type Array.Design.REF
## 164_I_.CEL biotin P-MTAB-41366 164_I_ array assay A-AFFY-141
## 164_II.CEL biotin P-MTAB-41366 164_II array assay A-AFFY-141
## 183_I.CEL biotin P-MTAB-41366 183_I array assay A-AFFY-141
## 183_II.CEL biotin P-MTAB-41366 183_II array assay A-AFFY-141
## 2114_I.CEL biotin P-MTAB-41366 2114_I array assay A-AFFY-141
## 2114_II.CEL biotin P-MTAB-41366 2114_II array assay A-AFFY-141
## Term.Source.REF Protocol.REF.4 Array.Data.File
## 164_I_.CEL ArrayExpress P-MTAB-41367 164_I_.CEL
## 164_II.CEL ArrayExpress P-MTAB-41367 164_II.CEL
## 183_I.CEL ArrayExpress P-MTAB-41367 183_I.CEL
## 183_II.CEL ArrayExpress P-MTAB-41367 183_II.CEL
## 2114_I.CEL ArrayExpress P-MTAB-41367 2114_I.CEL
## 2114_II.CEL ArrayExpress P-MTAB-41367 2114_II.CEL
## Comment..ArrayExpress.FTP.file.
## 164_I_.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 164_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 183_I.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 183_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 2114_I.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## 2114_II.CEL ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip
## Factor.Value.disease. Factor.Value.phenotype.
## 164_I_.CEL Crohn's disease non-inflamed colonic mucosa
## 164_II.CEL Crohn's disease inflamed colonic mucosa
## 183_I.CEL Crohn's disease non-inflamed colonic mucosa
## 183_II.CEL Crohn's disease inflamed colonic mucosa
## 2114_I.CEL Crohn's disease non-inflamed colonic mucosa
## 2114_II.CEL Crohn's disease inflamed colonic mucosa这么展示更利于查看
t(Biobase::pData(raw_data)[1,])
## 164_I_.CEL
## Source.Name "164_I"
## Characteristics.individual. "164"
## Characteristics.organism. "Homo sapiens"
## Characteristics.disease. "Crohn's disease"
## Characteristics.organism.part. "colon"
## Characteristics.phenotype. "non-inflamed colonic mucosa"
## Material.Type "organism part"
## Protocol.REF "P-MTAB-41361"
## Protocol.REF.1 "P-MTAB-41363"
## Extract.Name "164_I"
## Protocol.REF.2 "P-MTAB-41364"
## Labeled.Extract.Name "164_I:Biotin"
## Label "biotin "
## Protocol.REF.3 "P-MTAB-41366"
## Assay.Name "164_I_"
## Technology.Type "array assay"
## Array.Design.REF "A-AFFY-141"
## Term.Source.REF "ArrayExpress"
## Protocol.REF.4 "P-MTAB-41367"
## Array.Data.File "164_I_.CEL"
## Comment..ArrayExpress.FTP.file. "ftp://ftp.ebi.ac.uk/pub/databases/microarray/data/experiment/MTAB/E-MTAB-2967/E-MTAB-2967.raw.1.zip"
## Factor.Value.disease. "Crohn's disease"
## Factor.Value.phenotype. "non-inflamed colonic mucosa"筛选并保留关心的样本属性信息，个体信息 (Source.Name,Characteristics.individual.)，疾病信息 (Factor.Value.disease.),
发炎与否 (Factor.Value.phenotype.)。
Biobase::pData(raw_data) <- Biobase::pData(raw_data)[, c("Source.Name", "Characteristics.individual.", "Factor.Value.disease.", "Factor.Value.phenotype.")]原始数据质控
exprs(raw_data)可获取原始的表达信息，行代表芯片探针在芯片上的位置，列代表每个样品。
Biobase::exprs(raw_data)[1:5, 1:5]
## 164_I_.CEL 164_II.CEL 183_I.CEL 183_II.CEL 2114_I.CEL
## 1 4496 5310 4492 4511 2872
## 2 181 280 137 101 91
## 3 4556 5104 4379 4608 2972
## 4 167 217 99 79 82
## 5 89 110 69 58 47对表达数据取对数，并进行主成分分析。每个点代表一个样品，颜色代表是否发炎，形状代表疾病类型。从图中可以看出，在第一主成分两种疾病区分比较明显，在第二主成分疾病的区别也略大于是否发炎。而我们的关注点是发炎，而不是不同类型的疾病，需要在下游分析中考虑移除疾病的影响。
exp_raw <- log2(Biobase::exprs(raw_data))
PCA_raw <- prcomp(t(exp_raw), scale. = FALSE)
percentVar <- round(100*PCA_raw$sdev^2/sum(PCA_raw$sdev^2),1)
sd_ratio <- sqrt(percentVar[2] / percentVar[1])
dataGG <- data.frame(PC1 = PCA_raw$x[,1], PC2 = PCA_raw$x[,2],
Disease = pData(raw_data)$Factor.Value.disease.,
Phenotype = pData(raw_data)$Factor.Value.phenotype.,
Individual = pData(raw_data)$Characteristics.individual.)
ggplot(dataGG, aes(PC1, PC2)) +
geom_point(aes(shape = Disease, colour = Phenotype)) +
ggtitle("PCA plot of the log-transformed raw expression data") +
xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) +
ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) +
theme(plot.title = element_text(hjust = 0.5))+
coord_fixed(ratio = sd_ratio) +
scale_shape_manual(values = c(4,15)) +
scale_color_manual(values = c("darkorange2", "dodgerblue4"))另外，需要看下原始数据的表达分布。从图中看出，不同的芯片原始表达分布不均一，需要做合适的标准化处理。
# oligo::boxplot(raw_data, target = "core",
# main = "Boxplot of log2-intensitites for the raw data")
dataBoxplot <- data.frame(t(exp_raw), Sample=colnames(exp_raw),
Disease = pData(raw_data)$Factor.Value.disease.,
Phenotype = pData(raw_data)$Factor.Value.phenotype.,
Individual = pData(raw_data)$Characteristics.individual.)
dataBoxplot <- reshape2::melt(dataBoxplot, id.vars=c("Sample","Disease","Phenotype","Individual"))
ggplot(dataBoxplot, aes(x=Sample, y=value)) +
geom_boxplot(aes(color=Disease)) + theme_bw() +
theme(axis.text.x=element_text(angle=45, hjust=1, vjust=1))arrayQualityMetrics可以对芯片进行更多层面的检测，并生成一个检测报告。
arrayQualityMetrics(expressionset = raw_data,
outdir = "ysx",
force = TRUE, do.logtransform = TRUE,
intgroup = c("Factor.Value.disease.", "Factor.Value.phenotype."))背景校正
鉴于探针的荧光强度会受到非特异性杂交和光学检测系统噪音的影响，因此需要对信号强度进行背景校正，从而能更真实的反映和比较特异性杂交产生的信号，更好的反应基因表达情况。
芯片间标准化 (calibration)
不同的芯片杂交受很多因素影响，比如反转录效率，标记或杂交反应，芯片物理性质，试剂批次和实验室环境等，校正后，才可以在不同芯片之间可比。
Summarization
在Affymetrix平台，一个转录本通常用多个探针检测。对每个基因来说，需要根据所有探针的信息估算出其相对于总的转录本的表达值，也就是Summarization。随后就可以使用注释数据集获得基因的symbols或ENSEMBL ID。我们这个平台的注释信息存储在hugene10sttranscriptcluster.db包中。
head(ls("package:hugene10sttranscriptcluster.db"))
## [1] "hugene10sttranscriptcluster"
## [2] "hugene10sttranscriptcluster_dbconn"
## [3] "hugene10sttranscriptcluster_dbfile"
## [4] "hugene10sttranscriptcluster_dbInfo"
## [5] "hugene10sttranscriptcluster_dbschema"
## [6] "hugene10sttranscriptcluster.db"Affymetrix芯片”probesets”的变更
在采用芯片检测样品前，需要先提取RNA，加上一个荧光标签，并进行至少一轮PCR扩增以增加检测量。传统上这一步是根据mRNA3'端的poly-A序列设计引物完成反转。这样最终的扩增产物里面3’端的序列会更多些。而且如果起始RNA片段化了或降解了，5’端被检测到的几率更低了。因此设计探针时，更多的探针也是针对3’端设置的，这就是3' IVT Affymetrix芯片。它是基于探针集的，探针集的一组固定探针用于检测一个特定的基因或转录本。在一些情况下，如选择性Poly-A位点出现时，一个基因需要设计多个探针集检测。
而后来的Gene/ExonAffymetrix芯片是基于外显子的，探针集来源于外显子区。这时为了获得基因的表达，需要做两步summarization，第一步是从探针集获得外显子的表达，然后根据多个来源于同一个基因或转录本的外显子获取基因的表达。因此需要合适的基因注释包，如我们这用到的hugene10sttranscriptcluster.db。
Gene芯片相比于Exon芯片探针更少，只保留Exon芯片的一部分”good”探针。在Exon芯片上至少需要4个探针来代表一个Exon,而在Gene芯片上，很多代表一个基因的探针集只有3个或更少的探针。
如下图所示，用单一颜色代表一个探针集，一个基因的探针集用一个颜色集合表示，如所有黄色探针代表一个外显子，所有黄色、橙色、红色探针集都用于检测同一个基因。
左图中，每个外显子或探针集包含很多探针
(同样颜色的块出现多次)，因此做探针集或外显子水平的summarization是合适的。而对于gene芯片，每个探针集的探针数目很少，因此在探针集或外显子水平的summarization是不推荐的，虽然也可以通过包hugene10stprobeset.db来完成。
而且在Gene和Exon芯片上不再有错配探针 (mismatch
probes)。错配探针本来是用做背景校正排除非特异性杂交和背景噪音的，但被更好的计算方式取代。
一步完成背景校正、标准化、表达整合
oligo包提供了一个函数rma允许一步完成基于去卷积的背景校正、分位数校正标准化
(quantile normalization)和RMA (robust multichip average) 表达整合
(summarization)。
相对log表达 (RLE) 质控分析
我们先利用rma函数进行背景校正和summarization获得基因的表达量，而略过标准化。rma函数输出的结果已经做过对数处理。
palmieri_eset <- oligo::rma(raw_data, target = "core", normalize = FALSE)
## Background correcting
## Calculating Expression然后进行RLE (Relative Log Expression)分析，首先计算每个转录本在所有样品中的表达中位数，然后每个转录本的表达减去这个中位数，整理格式进行绘图。
row_medians_assayData <-
Biobase::rowMedians(as.matrix(Biobase::exprs(palmieri_eset)))
# 减去中值
RLE_data <- sweep(Biobase::exprs(palmieri_eset), 1, row_medians_assayData)
RLE_data <- as.data.frame(RLE_data)
# 转换为长格式
RLE_data_gathered <-
tidyr::gather(RLE_data, patient_array, log2_expression_deviation)
ggplot2::ggplot(RLE_data_gathered, aes(patient_array, log2_expression_deviation)) +
geom_boxplot(outlier.shape = NA) +
ylim(c(-2, 2)) +
theme(axis.text.x = element_text(colour = "aquamarine4",
angle = 60, size = 6.5, hjust = 1 ,
face = "bold"))图中，Y-轴代表每个芯片中转录本表达与该转录本在所有芯片表达中位数的偏差。箱体越延展表明芯片中转录本的表达量与中位表达差别越大，箱体在Y-轴的偏移表示大部分转录本的表达定量存在系统的增高或降低。这可能存在质量问题或是检测批次效应造成的。因此如果一个样品的boxplot的形状和中值与其它样品差别很大，需要进一步评估是否为异常样本，并考虑移除。
图中5个中位数在-0.7~-0.8之间的样品 2826_I, 2826_II, 3262_II,3302_II和3332_II可能是异常样品。需要结合后面的聚类分析，再确认是否移除。
校正三部曲
background-correct, normalize and summarize.
palmieri_eset_norm <- oligo::rma(raw_data, target = "core")
## Background correcting
## Normalizing
## Calculating Expression参数target定义summarization的程度，默认是core表示获得基因水平的表达
(using transcript clusters containing “safely” annotated genes)
(对gene芯片，只有这么一个选项)。如果是外显子芯片，还可以选择extended和full。如果想获得外显子水平的表达，也可以用probeset做参数，但是一般不推荐。
除了RMA还有其它的背景校正和标准化方法，但RMA通常是一个好的默认选择。RMA在所有芯片间使用quantile normalization标准化方法使得不同样品检测的表达值可比。
出校准化后的表达值
norm_exprs <- Biobase::exprs(palmieri_eset_norm)
norm_exprs <- data.frame(ID=rownames(norm_exprs), norm_exprs)
write.table(norm_exprs, "normalized_probe_expr.tsv", sep="\t", row.names=F, quote=F)校准化后数据的质量评估
主成分分析
exp_palmieri <- Biobase::exprs(palmieri_eset_norm)
PCA <- prcomp(t(exp_palmieri), scale = FALSE)
percentVar <- round(100*PCA$sdev^2/sum(PCA$sdev^2),1)
sd_ratio <- sqrt(percentVar[2] / percentVar[1])
dataGG <- data.frame(PC1 = PCA$x[,1], PC2 = PCA$x[,2],
Disease = Biobase::pData(palmieri_eset_norm)$Factor.Value.disease.,
Phenotype = Biobase::pData(palmieri_eset_norm)$Factor.Value.phenotype.)
ggplot(dataGG, aes(PC1, PC2)) +
geom_point(aes(shape = Disease, colour = Phenotype)) +
ggtitle("PCA plot of the calibrated, summarized data") +
xlab(paste0("PC1, VarExp: ", percentVar[1], "%")) +
ylab(paste0("PC2, VarExp: ", percentVar[2], "%")) +
theme(plot.title = element_text(hjust = 0.5)) +
coord_fixed(ratio = sd_ratio) +
scale_shape_manual(values = c(4,15)) +
scale_color_manual(values = c("darkorange2", "dodgerblue4"))第一主坐标轴大体可以分开发炎与未发炎样品，第二主坐标轴分开两种不同的疾病。(从图上看，还是疾病的区分更明显)
热图聚类分析
通过热图展示样品之间的相似度关系，并且标记样品的属性如是否发炎或是否来自同一个疾病。
# 生成行注释
phenotype_names <- ifelse(str_detect(pData(palmieri_eset_norm)$Factor.Value.phenotype., "non"), "non_infl.", "infl.")
disease_names <- ifelse(str_detect(pData(palmieri_eset_norm)$Factor.Value.disease., "Crohn"), "CD", "UC")
annotation_for_heatmap <- data.frame(Phenotype = phenotype_names, Disease = disease_names)
row.names(annotation_for_heatmap) <- row.names(pData(palmieri_eset_norm))计算Manhattan距离，并进行聚类分析
# 使用dist计算样本之间的距离，首选做下转置，因为dist是对列进行操作的
# dist默认是欧式距离，反应直线路径的平方距离
# 这里用Manhattan距离，计算的是直角路径的绝对距离，计算结果更稳定
dists <- as.matrix(dist(t(exp_palmieri), method = "manhattan"))
rownames(dists) <- row.names(pData(palmieri_eset_norm))
hmcol <- rev(colorRampPalette(RColorBrewer::brewer.pal(9, "YlOrRd"))(255))
colnames(dists) <- NULL
# 设置对角距离为NA，增强颜色的展示力
diag(dists) <- NA
# www.ehbio.com/ImageGP
ann_colors <- list(
Phenotype = c(non_infl. = "chartreuse4", infl. = "burlywood3"),
Disease = c(CD = "blue4", UC = "cadetblue2"))
pheatmap(dists, col = (hmcol),
annotation_row = annotation_for_heatmap,
annotation_colors = ann_colors,
legend = TRUE,
treeheight_row = 0,
legend_breaks = c(min(dists, na.rm = TRUE),
max(dists, na.rm = TRUE)),
legend_labels = (c("small distance", "large distance")),
main = "Clustering heatmap for the calibrated samples")热图中上部黄色的条带可能是异常样品 (2826_II, 3262_II, 3271_I, 2978_II 和 3332_II)，部分样品与之前RLE图推测可能是异常的样品重合，可以考虑移除。这里为了跟原文一致，所有样品都保留了下来。
另外arrayQualityMetrics有更精确的计算异常样品的方式，我们后续也会提供基于网络图距离的异常值剔除方法。
过滤低表达基因
芯片数据中有一大部分探针的信号值与背景值相差不大，而且在不同样品中差别不大，不能给我们的生物问题提供更多信息。(注：虽然有部分表达低的基因表达变化比较小时就可以有比较大的作用，但因为其表达低，检测噪音也大，结果可信度一般)
首选绘制下每个转录本在所有芯片中的表达中位数的分布
(共计33,297转录本)。可以看到左侧区域富集低表达的转录本，具体选择移除的阈值主观性比较强，这里选择表达中位数需要大于4，给后续分析多保留一些基因。也可以排序一下，保留前75%的基因。
palmieri_medians <- rowMedians(Biobase::exprs(palmieri_eset_norm))
hist_res <- hist(palmieri_medians, 100, col = "cornsilk1", freq = FALSE,
main = "Histogram of the median intensities",
border = "antiquewhite4",
xlab = "Median intensities")man_threshold <- 4
hist_res <- hist(palmieri_medians, 100, col = "cornsilk", freq = FALSE,
main = "Histogram of the median intensities",
border = "antiquewhite4",
xlab = "Median intensities")
abline(v = man_threshold, col = "coral4", lwd = 2)假如最小的实验组有5个样品，如果一个转录本检测到表达值大于之前设定的阈值
(man_threshold)的样本数小于5，则移除。
首先获得各个实验组的样品数
no_of_samples <-
table(paste0(pData(palmieri_eset_norm)$Factor.Value.disease., "_",
pData(palmieri_eset_norm)$Factor.Value.phenotype.))
no_of_samples
##
## Crohn's disease_inflamed colonic mucosa
## 15
## Crohn's disease_non-inflamed colonic mucosa
## 15
## ulcerative colitis_inflamed colonic mucosa
## 14
## ulcerative colitis_non-inflamed colonic mucosa
## 14
# 获得最小样品数
samples_cutoff <- min(no_of_samples)
# sum(x > man_threshold)：转录本在多少样品中表达值高于给定阈值
idx_man_threshold <- apply(Biobase::exprs(palmieri_eset_norm), 1,
function(x){sum(x > man_threshold) >= samples_cutoff})
table(idx_man_threshold)
## idx_man_threshold
## FALSE TRUE
## 10493 22804
# 提取过滤转录本
palmieri_manfiltered <- subset(palmieri_eset_norm, idx_man_threshold)注释转录本簇
在我们构建线性模型做差异分析之前，先把探针映射到Gene symbol，并只保留这部分基因用于后续分析。
# 使用select函数，获得探针ID对应的Gene symbol和描述
anno_palmieri <- AnnotationDbi::select(hugene10sttranscriptcluster.db,
keys = (featureNames(palmieri_manfiltered)),
columns = c("SYMBOL", "GENENAME"),
keytype = "PROBEID")
# 过滤掉没有Symbol的探针
anno_palmieri <- subset(anno_palmieri, !is.na(SYMBOL))
head(anno_palmieri)如果一个探针对应于多个Gene symbol，我们没有办法判断到底是哪个，只能移除。
# 按探针ID分组
anno_grouped <- group_by(anno_palmieri, PROBEID)
# 计算每组Symbol的个数
anno_summarized <-
dplyr::summarize(anno_grouped, no_of_matches = n_distinct(SYMBOL))
head(anno_summarized)
# 过滤Symbol个数大于1的
anno_filtered <- filter(anno_summarized, no_of_matches > 1)
head(anno_filtered)
probe_stats <- anno_filtered
nrow(probe_stats)
## [1] 1763有1763个探针ID对应多个Gene Symbol，从Assya数据中移除这些探针ID。
# featureNames(palmieri_manfiltered): 获取所有探针名字
ids_to_exlude <- (featureNames(palmieri_manfiltered) %in% probe_stats$PROBEID)
table(ids_to_exlude)
## ids_to_exlude
## FALSE TRUE
## 21041 1763
palmieri_final <- subset(palmieri_manfiltered, !ids_to_exlude)
validObject(palmieri_final)
## [1] TRUE前面是从assay数据中移除了这些探针ID，后面需要从feature data中移除:
head(anno_palmieri)现在的feature data是空的，增加一列PROBID。
# 增加一列PROBID
fData(palmieri_final)$PROBEID <- rownames(fData(palmieri_final))
# 使用left_join根据PROBID列，合并，增加注释信息
fData(palmieri_final) <- left_join(fData(palmieri_final), anno_palmieri)
## Joining, by = "PROBEID"
# restore rownames after left_join
rownames(fData(palmieri_final)) <- fData(palmieri_final)$PROBEID
validObject(palmieri_final)
## [1] TRUER统计和作图
随机森林randomForest 分类Classification 回归Regression
GEO是当今最大、最全的公共基因数据资源库，包括基因的表达、突变、修饰等信息，涵盖几乎所有的疾病，且单个实验检测样品数目较多。TCGA数据库包含11,000个病人的33种肿瘤的7个不同层面的基因数据 (包括基因表达、CNV，SNP，DNA甲基化，miRNA，外显子组等)和临床数据，意在解析癌症发生的分子机制、肿瘤的亚型和治疗靶点等。
这两个来源的数据都是对外公开的，是学习、研究和应用的一个资源宝库。从2006年TCGA计划启动以来，基于TCGA数据发表的文章呈指数增加，一大部分来源于对TCGA数据的再次挖掘。因此学习利用生物信息技术挖掘GEO/TCGA公共数据中疾病的分子特征、合适的检测指标具有重要的临床和科研价值。本课程将从GEO/TCGA的表达、突变数据入手，探索公共数据挖掘的基本套路，分享数据分析和可视化的思路和代码，以便应用于自己的研究。
课程涉及主要内容如下：
每节课1小时一个主题，理论结合实战，学懂原理，实战实操，全是老司机多年经验和代码的无私分享。利用自己电脑，轻松实现整套分析。如果有基础，可以多理解代码内容，做更多定制。如果基础弱一些，只需修改几个备注的变量，即可完成全部分析。下面是课程安排，如11代表第一天第一节课，26代表第二天第六节课，(实际上课会有调整，理论和实战穿插调和)，41为两周后的线上集中视频答疑。
| 编号 | 主题 | 简介 | |
|---|---|---|---|
| 11 | GEO挖掘案例介绍和结果解读 | 学习套路和宏观把控 | |
| 12 | R语言基础知识介绍 | 基础变量和数据操作 | |
| 13 | ggplot2绘图基础 | 热图、火山图等常见图绘制 | |
| 14 | GEO数据库使用介绍 | 数据查找、下载、清洗、批次校正、可视化 | |
| 15 | 芯片全套分析 | 差异基因、GO GSEA(时间序列)富集分析、可视化 | |
| 16 | WGCNA共表达网络 | KEGG/Reactome通路可视化 | |
| 21 | 实战重现2018 纯公共数据Science文章 | 文章整体解读和亮点分析 | |
| 22 | 实战重现2018 纯公共数据Science文章 | 表达模式评估，样品差异分析 | |
| 23 | 实战重现2018 纯公共数据Science文章 | 不同来源数据校正和比较的原理和操作 | |
| 24 | 实战重现2018 纯公共数据Science文章 | WGCNA共表达模块分析、网络可视化、模块性状关联 | |
| 25 | 实战重现2018 纯公共数据Science文章 | 基因表达和突变数据关联 | |
| 26 | 常见图快速绘制、解读和排版 | 见图 | |
| 31 | TCGA数据 | 案例介绍和结果解读 | |
| 32 | 实战重现2018 JAMA文章 | TCGA数据下载和整理 | |
| 33 | 实战重现2018 JAMA文章 | TCGA数据表达分析全套 | |
| 34 | 实战重现2018 JAMA文章 | 突变模式, 突变负荷和生存分析 | |
| 35 | 实战重现2018 JAMA文章 | 突变与临床因素相关性分析以森林图展示 | |
| 36 | 考试、圆桌论坛 | 自评学习效果、知识点回顾 | |
| 41 | 答疑-线上 | 答疑、考试内容串讲 |
R基础知识和图形绘制
GEO数据分析
GEO数据库使用介绍, 数据查找、下载、清洗、批次校正、可视化
芯片全套分析, 差异基因、GO GSEA(时间序列)富集分析、可视化
实战重现纯公共数据Science文章, 整体解读, 亮点分析, 结果重现
TCGA数据分析
TCGA数据，案例介绍和结果解读
TCGA数据下载和整理，临床数据获取和整理
TCGA数据表达分析全套，差异基因富集分析等
实战重现2018 JAMA文章, 突变模式, 突变负荷和生存分析, 突变与临床因素相关性分析以森林图展示, 新队列数据验证结果
希望大家报名时，给出自己想重复的文章或结果，我们综合评估，优先照顾，定制专属课程。如果当时不能满足的，也会在后续讨论群提供解决方案，毕竟我们为所有线下学员提供免费绘图。
(如果基础薄弱，报名付款成功后，可免费领取基础程序课，做好准备工作, 让程序成为我们的得力工具而不是学习新知识的绊脚石。)
往期精彩回顾
主讲教师
陈同，博士，2015毕业于中科院遗传与发育生物学研究所，生物信息专业博士，在Cell Stem Cell(IF=23.2，第一作者兼封面文章)，Nucleic Acids Research，Stem Cells and Development等高水平杂志以第一作者或主要作者发表文章，运营有数万人关注的《生信宝典》微信公众号，给你不一样的学习生信体验。
刘永鑫，博士。2008年毕业于东北农大微生物学专业。2014年中科院遗传发育所获生物信息学博士学位，2016年博士后出站留所工作，任宏基因组学实验室工程师，目前主要研究方向为宏基因组学数据分析与可重复计算。发表论文10余篇，SCI收录7篇。2017年7月创办“宏基因组”公众号，目前分享宏基因组、扩增子原创文章185篇，关注人数2.3万人，累计阅读近300万次。
授课模式
本课程以讲解流程和实际操作为主，采用独创四段式教学，封装好的代码全部分享，随处可用：
第一阶段 3天集中授课；
第二阶段 自行练习2周；
第三阶段 在线直播答疑；
第四阶段 培训视频继续学习；
实现教-练-答-用四个环节的统一协调。
培训时间
2019-04-19 到 2019-04-21 (线下讲解实战)
每天早9点到晚6点，半封闭式教学 (最后1小时为集中讨论时间，最后一天会稍微提前一些，多留出时间讨论，也方便老师乘车返回)
报到时间：培训当天
授课地点
北京市西城区鼓楼附近 (具体开课前一周通知)
课程价格
截止 2019-04-12前 4500 元/人
名额有限，每次课程报名满40人后自动关闭报名通道
提供易汉博基因科技实习机会或工作机会
课程连报有优惠，具体见报名网站
课程福利
座位按报名并缴费或预付款成功顺序从前到后龙摆尾式排序
赠送程序基础课一份 (http://bioinfo.ke.qq.com)
多人 (N，10>N>1) 组团报名并同时缴费，每人还可减免N-1百元 (最高500)
赠送金士顿U盘一个（32G含培训数据和脚本）
附推荐语分享对应的招生信息到朋友圈，截图发到train@ehbio.com 可获得200元生信宝典腾讯课堂课程优惠券（可拆分供多个课程使用）
注意事项 *
需自备笔记本电脑，推荐使用win10系统，4G以上内存 (推荐8G)。课程实践根据需要会提供云计算平台
培训班所有数据，文档为内部资料，仅供参阅，未经允许不得翻印外传登刊
上课期间禁止录音，录像
成功付款的学员，若临时有紧急事情不能到来的，可申请延期，更换后续培训班；也可申请退款
若开课2周 (含) 前申请退款可退还85%费用；开课3个工作日 (含) 前申请退款退还70%的费用 (若已开发票需承担相应手续费)
不可先延期再退款
更多课程的详细介绍，请扫描下方二维码。
复制以下链接http://www.ehbio.com/Training/ 或 点击阅读原文跳转报名页，成为实验中不可或缺的人，赶快报名吧！
易生信系列培训课程，扫码获取免费资料
更多阅读
后台回复“生信宝典福利第一波”获取教程合集
听说分享到朋友圈的朋友会在公众号周年庆时中奖 (大家还记得去年的大放送吧，不记得查查历史